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Abstract: Forthcoming large angle surveys are planned to obtain high precision tomo- 
graphic shear data. In principle, they will allow us to recover the spectra of matter density 
fluctuation, at various redshift, through the inversion of the expressions yielding shear from 
fluctuation spectra. This was discussed in previous work, where SVD techniques for matrix 
inversion were also shown to be the optimal tool to this aim. Here we show the significant 
improvements obtainable by using a 7 bin tomography, as allowed by future Euclid data, 
as well as the question of error propagation from shear to fluctuation spectra. We find that 
the technique is a promising tool, namely for the analysis of baryon physics through high-£ 
shear spectra and to test the consistency between expansion rate and fluctuation growth. 



hydrodynamic simulations. 
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1. Introduction 



Dark Energy (DE) is one of the main discoveries - and puzzles - in contemporary physics. 
Cosmic Microwave Background spectra or Baryonic Acoustic Oscillations directly constrain 
its contribution to the cosmic budget, but only simultaneous low-z measures of cosmic 
expansion rate and fluctuation growth can provide us real clues on DE nature |TJ. In fact, 
besides of its state equation w(z), we aim at knowing whether DE is a further physical 
dark component, or its observations are a signal of new physics, as GR violations El, 
interactions between the dark components, possibly suggesting a unified picture for them 
H |4|, or even more exotic scenarios [||, |j. It is then worth outlining that, while SNIa data 
constrain the expansion rate, tomographic shear data can enable us to follow the growth of 
density inhomegeneities, matching it with background geometry, namely if such data can 
be suitably shared in a number of redshift bins. Moreover, being directly sensitive to the 
whole mass distribution, shear data allow us to forget any problem related to light-mass 
conversion. 
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In order to exploit this kind of data, it has become customary to follow a Bayesian 
approach. We then consider a set of parameters spanning a wide model set, for each 
of whom a number of observables can be predicted; prediction are then compared with 
observational data and relative errorbars, so obtaining likelihood hyper-elipso'ids in the 
parameter space. Similarly, when cosmic shear data will be available, one shall predict 
the tomographic shear spectra Cy(^) for the models considered (an operational definition 
of Cij{€) is given below; let us soon outline, however, that the indeces i,j run on the 
tomographic bins) and enrich the fit to data by comparing them with tomographic shear 
data. 

In a previous paper Paper I herebelow) we considered a somehow different option: 
directly deriving the fluctuation spectra P(k, z) at various redshift z from the tomographic 
shear spectra Cij{£). If doing so, a prediction of shear spectra is no longer required, and 
shear data can be soon compared with other data on density fluctuation growth, to be then 
used in parallel with them. This option makes sense, in particular, if we wish to inspect 
background geometry and fluctuation dynamics separately. 

In the literature, an integral relation yielding C\j(£) from P(k, z) has been known since 
long. A procedure to invert it has been introduced in Paper I. There, however, we restricted 
our analysis to 5 tomographic bins, while Euclid 1 data, obtained from a celestial area of 
15,000 square degrees with a median redshift z m = 0.9, will allow N-bin tomography, with 
N> 5 and up to 10 §§. 

In this paper, we consider a 7 bin tomography and find that the improvements allowed 
by this N increase, when aiming to recover P(k,z), are substantial. The recovery still uses 
the singular value decomposition (SVD) technique for matrix inversion. Specific options of 
this technique, suitable to treat quasi-singular matrices, will be however used here for the 
first time in this context. Furthermore, we shall treat the question of noise propagation 
from shear to fluctuation spectra. Matrix theory allows us to predict upper limits to error 
magnification in such transition. Lukily enough, in this specific case, such limits are only 
marginally approached, while a simple filter allows us further noise reduction. 

Let us also notice that, besides of verifying the coherence between expansion and 
fluctuation growths, this inversion technique could also facilitate the discrimination between 
different options on baryon physics, shaping P(k, z) at high-Zc. 

In order to test the inversion algorithm we need to input fluctuation spectra for a given 
model. We use them to build fluctuation spectra P(k, z) at various redshifts and then, from 
these spectra, the tomographic shear spectra Cij{£). Successively, starting from dj(£), 
we test how efficiently P(k, z) is recovered, under various options. Being also interested 
in baryon physics and, therefore, having to explore the high-£ region, hydrodynamical 
simulations are to be used. A brief discussion on the simulations used will be therefore 
provided. Let us however preliminarly outline that the lensing spectra can be roughly 
shared in 3 ^-ranges (see Figure |3|, here below): The C- L j{£) for £ <~ 500 essentially feel 
just the linear dynamics, even for i, j = 1. For 500 < £ < 1500 the contribution coming from 
non-linear k 's becomes relevant. At £ > 2000, shear data start to exhibit a dependence on 
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baryon physics. 

The plan of the paper is as follows: In the next Section we shall review the procedure 
allowing to pass from data on galaxy ellipticity to shear spectra. In particular, we shall 
obtain the window function in the 7-bin case, taking into account that, for most lensed 
galaxies, only the photometric redshift can be given. 

In Section 3 the fluctuation and shear spectra will be defined and their relation set in 
an operational form. Section 4 consists of three discussions: on the simulation used; on the 
pattern followed to interpolate its outputs, so obtaining P(k, z) at the prescribed k and z 
values; on the procedure of Gauss-Laguerre integration to pass from fluctuation to shear 
spectra. Section 5 introduces the formal procedure to invert the equation yielding Cij(£) 
from P(k,z), also based on the Gauss-Laguerre procedure. Section 6 then introduces 
the SVD technique and uses it for the study of the singularity level, as well to recover 
P(k, z) from Cij(£), in 3 cases: (i) When Cij(£) is formally obtained with a GL integration 
technique, (ii) When Cij(£) is "exact", (iii) When noise is added to it. Results are 
summarized and discussed in Section 7 and conclusions are drawn in Section 8, where we 
also outline possible options to improve the present approach. 



2. Evaluation of shear— shear correlation functions 



The procedure to derive shear spectra starts from approximating galaxy images by ellipses 
of axes a and b (a > b), and setting 

7 = 7+ + *7x ee 4tS c ^' 



in agreement with ||, [10|. Here ip is the angle between the major axis of the ellypse and 



the x-axis of a chosen coordinate system on the sky, assumed to be locally flat. When the 
galaxy sample occupy a wide celestial area, as in the case of Euclid, it can be necessary 



to correct for the sky curvature [11|. 



Measuring galaxy ellipticity is a non-trivial task, but the main problem is that ellip- 
ticities are not due to gravitational shear only. The so-called intrinsic shear due to nearby 
galaxy alignment is the main problem while, aiming at correlating galaxy ellipticities, the 
presence of an uncorrelated random shear component is not so critical. Assuming that 
data can be simply and reliably cleansed from intrinsic alignement may be premature (see 



however [12], [13] and references therein), but here we let apart this problem and proceed 
as though 7 were due to gravitational shear only. 

Measures will be considered in the context of (spatially flat) models whose background 
metric reads 

ds 2 = a 2 (r)[dr - d\ 2 } , (2.2) 

so that r is the conformal time, dX being the comoving distance element, and a(r) 
i 1 + z)~ l the scale factor. If tq is the present conformal time, 



U(z) = T - t(z) 



(2.3) 



is the conformal time distance from the redshift z. Owing to eq. fl2.2j ) it is also the comoving 
distance from z. 

According to data, we then assume that galaxies observed in unit solid angle, on the 
light cone of this space-time, have a normalized distribution on redshift reading 



n(z) 



d 2 N 
dfl dz 



exp 



with 
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(2.4) 



(2.5) 



and A = 2, B = 1.5, so that C = 1.5/zq (with zq = z m /1.412 obtained from the median 
redshift z m = 0.9 , in agreement with Euclid specifications ||). These phaenomenological 
values can however be suitably modified, if needed. 

In order to appreciate the effects of fluctuation evolution, observed galaxies are then 
shared into iV redshift bins, with limits z r possibly selected so that each of them contains 
an equal number of galaxies. However, if we want the galaxy set to be large, we must 
accept that only their photometric redshifts are available. 

To evaluate the expected distribution on redshift for the r-th bin galaxies, we must 
then apply the filter 
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(2.6) 



V2a(z) J ^ [\/2a(z)_ 

to n(z) (through this paper we shall take cr(z) = 0.05 (1 + z), coherently with EUCLID 
expectations pL see also so obtaining the distribution 



D r (z) = n(z)U r (z) . 



(2.7) 



whose integral is ~ 1/N, while their sum is (exactly) n(z), as shown in Figure |T[ 

The images of the galaxies belonging to a bin are lensed by the gravity of matter laying 
at lower z. This is taken into account by defining the window functions 



W r (z) = -O m (l + z) dz' S r (z') V 



Az r 



u(z') — u(z) 
ujz~>) 



here u(z) = tq — t(z) and 



S r (z) 



DJz) 



(2.8) 



(2.9) 



J °° D r (z')dz' 

is obtained by normalizing each D r to unity, and the integration can be restricted to the 
interval Az r , where D r is significantly > 0; let also be V{x) = x or if x > or < 0: only 
systems closer than a galaxy can distort its image. Window functions for the 7-bin case 
are shown in Figure 0. 
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2000 4000 6000 

( T o- T )/Mpc 

Figure 2: Window functions for the 7-bin case. In the spatially flat models considered, along the 
light cone, conformal time delay and comoving distance coincide. 



In order to calculate the window functions W r though the expression ( |2.8| ) we need to 
know the u{z) dependence, i.e., the background geometry. More in general, the reference 
cosmology in this paper is a ACDM model with density parameters for total matter and 
baryons, Hubble parameter, primordial spectral index and the mean square amplitude of 
linear density fluctuations, at z = 0, on the scale of 8 /t~ 1 Mpc, shown herebelow in Table I: 
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Tabic I 



il m Vt b h n s a s 

0.24 4.13 x 1(T 2 0.73 0.96 0.8 

However, as far as building W r is concerned: (i) no knowledge of baryon physics is needed, 
not even the value; (ii) no parameters ruling the dynamics, as as or n s , are to be known. 
The cosmology described by Table I is the same of Paper I. 

The shear-shear correlation functions among galaxy pairs can then be estimated through 
the sums 



7+;a7+ ; /3 +/ 7x;a7 X ;/3 



e 



a/3 



e°,_(0) = — ^ ^ (2.10) 

(a and f3 run over the galaxy set), by sharing the angle 9 into bins, and setting = 1/0 
when the angle 9 a p between galaxy images is within/outside the suitable bin. 

3. Density fluctuation and shear spectra 

The shear-shear correlation functions can be expanded in spherical harmonics on the ce- 
lestial sphere. It can then be shown that 

roc 

eL_{e) = (27T)- 1 / di£j 0/4 (ee)c tJ (£) , (3.1) 

Jo 

with the same angular spectrum Cij{£) {J n being first kind Bessel functions of order n). 
Owing to the finiteness of the celestial sphere, the modulus of a 2-dimensional wave vector £ 
would be allowed just discrete values. However, the contribution from the ^-th component 
is roughly due to shear on the angular aperture -d ~ 27r/^, and we do not expect shear 
signals to be significant at large-??. 

Here, as usual, the Fourier transform of the density fluctuation field e(x, z) = [p(x, z) — 
p(z)]/p(z) reads 

8(k,z) = — ^ / tfW k - x e(x,z) (3.2) 



(2vr) 3 d 

and, at any redshift z, we define the power spectrum P(k,z) = (\5(k, z)\ 2 ) . Then, it can 
be shown that the relation 



"TO 
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Cij(£) = H 4 / du Wi(u)Wj(u) P[£/u,u] (3.3) 
J 

holds. Here and in what follows, both P(k, z) and P(k, u) will be used, the latter symbol 
indicating P[£/u, z(u)}, as in eq. (p.3|). The Wi(u) are the window functions defined in 
eq. ( |2.8| ). Here again notice the rough correspondence k ~ 2ir/\. 

It is reasonable to expect that the contributions to the shear on an angular scale $ 
depend on density fluctuations over the increasingly large scale A subtended by as z 
increases. Accordingly, eq. (^^) shows that the shear signal at a given £ receives contri- 
butions from decreasing k values as u increases. In principle, for u — > 0, P(k, u) should be 
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Figure 3: The red lines show the k values contributing to the Cij(£) shear spectra at constant £ 
(aside them) and growing redshift. 

evaluated at k — > oo where, however, it vanishes. 

If we consider the log fc-log u plane, the integration in eq. ( [3 .3D is carried along tilted 
straight lines as those shown in Figure ^, each line corresponding to a given £. 

The spectra P(k,z) therefore yield functions Pe(u) = P[£/u, z(u)] with £ fixing a 
straight line on the log /c-log u plane and u being an abscissa on such line. 

4. Prom fluctuation to shear spectra 

Our aim now amounts at finding a technique allowing the inversion of eq. (|3.3| ). To do 
so, we shall first work out the spectra C^{€) for a given model, according to eq. ( |Q| ) 
itself. In the linear range there is no problem to obtain the function Pe(u) = P[£/u, z(u)], 
for any u and £, from algorithms like CAMB. On the non-linear range one has then to 
make recourse to simulations, which however yield P(k = £/u,u), on a set of points kb~u a 
(redshift values z a ). 

4.1 Spectra interpolation 

There is then an interpolation problem, to work out Pe(u) in the points needed for inte- 
gration. While integration routines are very efficient and allow any required precision, the 
actual precision depends on how interpolation is performed among the limited set of points 
available. 

To do it, there are at least 3 different patterns, schematically shown in Figure ||. We 
can first work out Pi{u a ) at fixed £, by interpolating among kf, values. To perform the 
integral we need Pi(u) for suitable values of u ^ u a , and this is obtainable by further 
interpolating at constant £. In Figure || the pattern 1 outlines this procedure. (It may be 
worth outlining that, here and in the next cases, interpolation does not only involve the 2 



Log(k) 



Figure 4: Interpolation patterns. This schematic Figure assumes that "data" are available at the 
crossings between horizontal and vertical dotted lines. The tilted lines are at constant t. Along one 
of them, we schematically indicate a set of points where we need to know P(l 7 u), to perform the 
numerical integration. To reach them, in any case, there are 2 interpolation steps. The direction 
of former (latter) is indicated by the two parallel segments (a single segment, tilted or orthogonal). 

closest points in each direction. The thick segments in the Figure are to be taken just as 
an indication of directions along which one uses a cubic spline.) 

Otherwise, we can first determine P(kb, u) at the u value needed and then we interpo- 
late at constant u along k (pattern 2). 

The third option considered is determining first P(k, u a ) at the k value needed; then 
we interpolate at constant k along u (pattern 3). 

In order to test which procedure is best, we made use of the HALOFIT algorithm, 
yielding spectra at arbitrary k and u. This allowed us to appreciate that the procedures 2 
& 3 are substantially equivalent. On the contrary, the procedure 1, requiring interpolation 
along the tilted ^-constant curves, yields consistent results only when spectra are known 
at many more redshift values z a (approximately 3 times more). 

A similar procedure was also used to test that the k and z number where we have 
simulation "data" is adequate; as a matter of fact, it turns out to be so only if we avoid 
the first interpolation pattern. 

The integration interval in eq. ( |3.3[) is apparently finite. When u reaches the conformal 
age of the Universe To, however, z approaches oo. Figure || shows that all W\ vanish well 
before so, and this sets an effective upper limit to the integration. 
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Figure 5: The magenta dashed line sets the k upper limit above which shot noise covers physical 
spectra. The black curves are the integration patterns for the I values indicated aside them (at vari- 
ance from the previous Figure, due to the different ordinate, here they are not straight lines). The 
horizontal dashed lines indicate are at z = 0.06, 0.5, 1. Accordingly, even the power spectra needed 
to obtain shear spectra at £ = 10000 are free of shot noise as soon as z >~ 0.06, corresponding to 
- 17.5 h^Mpc. 

Aiming at 5— digit precision, a numerical integration performed by summing on a large 
number of equispaced points requires ~ 2000 points up to u ~ 6000 (z ~ 3), paying atten- 
tion that each contributions added to the sum is sufficiently large in respect to numerical 
precision. 

4.2 The simulation 



Let us give now some information on the simulation used. More details are given in [15]. 
The simulation was run in a box of side L = 410 Mpc (corresponding to k^ ~ 1.5 x 
10" 2 /iMpc -1 ) where (2x)1024 3 particles were set, and with a force resolution e = 7.5 h 
kpc (corresponding to k € ~ 8.4 x 10 2 /iMpc _1 ). More precisely, spectra can be used up 
to k ~ N(2tt/L) with N ~ 2 15 = 32768, however keeping k (~ 5.0 x 10 2 /iMpc" 1 ) < k e , 
provided that they are not covered by shot noise (see Figure ||). 

Being hydrodynamical, two populations of particles are used in simulations. With 
the cosmological parameter of Table I, the masses of CDM and baryonic particles are 
m c ~ 1.89 x 10 9 h~ l M & and mt, ~ 3.93 x 10 8 h~ l M & , respectively. The baryon dynamics 
includes the effect of cooling and star formation. A description of metal production from 



chemical enrichment contributed by SN-II, SN-Ia and AGB stars, as described in [16], is also 
included. Stars of different mass, distributed according to a Salpeter IMF, release metals 
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over the time-scale determined by the corresponding mass-dependent life-times. Kinetic 
feedback is also implemented by mimicking galactic ejecta powered by SN explosions. In 
these simulations AGN feedback is not included, at variance from []l7]], where smaller 
simulation boxes are however kept. Simulation start at z = 41 and spectra are obtained 
for the redshift values 

l + z r = W r/20 (r = 1, ....,19) . (4.1) 

8 of them, up to z ~ 1.24 are between and 1.4; 7 more are within z = 4.02. These 15 
spectra are used for interpolation. 

4.3 Gauss Laguerre integration procedure 

Besides of performing "exact" Riemann integration, here we need to consider a Gaussian 
integration procedure as well. It amounts to projecting each integrand function f(u) = 
Wi(u)Wj(u) Pi(u) onto polynomials 7T a (u), orthogonal with an assigned weight function 
R(u), finding its components f a . As the integrals II Q of each polynomial are known, the 
integral of f(u) reads Yl a =i fcJ^ai an d this is a reliable value if N is large enough. For 
practical aims, this technique can be translated into a simple procedure, i.e. to a weighted 
sum of values taken by the integrand function f(u) in a suitable set of points u a . 
More in detail, using monic polynomials, we have that 

dx R{x) 7r ct (x)7r / 3(x) =M5 a p , (4.2) 





with a known normalization J\f. Monic polynomials are obtained from a suitable recurrence 
relation assuming that the coefficient of the leading term, for each a, is unity. If we 
then truncate the sum to N terms, the N zero's of itn(x) are the points x a , while the 
corresponding weights are 

/ °°rfs i?(a) 7 r^_ 1 (x) 
w a = - ( — r-rt — s— , (4.3) 

n'(x) being the ordinary derivative of n(x). It is then 

\ duf(u)= dxF{x)R{x) = y2w a F(x a ) (4.4) 
Jo Jo a 

with 

F(x) = f[u{x)](du/dx)/R(x) . (4.5) 

For R(x) oc e~ x , this technique is dubbed Gauss-Laguerre integration, as TT a (x) = L a (x), 
the Laguerre polynomials. In the case of eq. ( |3.3| ), where integration is cut off by the decay 
of the Wi functions, this approach can be applied by assuming x = (u/u) 13 and suitably 
selecting then u and j3. 

5. Formal inversion 

Let us then rewrite eq. ( |3.3j ) as follows: 

N N 



Here we have set 



A = ij, c A = Cij/Ho 



(5.2) 



with the correspondence law 



hJ 
A 



1,1 



1 



1,7 2,2 
7 8 



2,7 3,3 
13 14 



7,7 
28 



while, by setting u r 



u(x r ), we have 



Wi (x r ) Wj (x r ) u r I [/3x r R ( x r ) ] 



(5.3) 



Pr{£) = Px r (Z) = PeM = P(t/u r ,u r ) . 



(5.4) 



If, in eq. (|5.1| ), we take N = 28, M.Ar are square matrices and, provided that they are not 
singular, the inverse equation 



also holds. In principle, we should then be able to work out the spectrum P(k, z) for any 
k = £/u r , at the redshift values z r = z(u r ). 

Notice that the inversion procedure acts on each £ value separately. Any z and/or k 
value can be attained, in principle, just by suitably choosing the u and f3 parameter. In 
principle, we could perform a different choice for each I. 

This procedure was first introduced in Paper I, where N could reach atmost 15. This 
led to 2 difficulties: (i) 15 gaussian points were far from allow a nearly-exact integration; 
(ii) the determinant of the 15 x 15 A4 matrix was close to singular. 

Curiously enough, an inversion could be attained by passing from 15 to 12 integration 
point. Although worsening the (i) point, this allowed to formulate a redundant system of 
equations, which could then be treated by using the SVD technique. A suitable choice of 
x and /3 allowed then to overcome the (ii) points and to achieve an exact inversion of the 
Cij(£) spectra formally obtained through Gaussian integration. 

As a matter of fact, however, even this was not yet enough when trying to invert the 
results of exact integration (i.e., the Riemann integration): some P(k, z) were recovered, in 
fact, but systematically ill-normalized. As a matter of fact, however, once the cosmological 
parameters needed to build the W r (u) are known, one needs to add just J)& and <7g to predict 
the exact normalization (the dependence on n s was found to be negligible). In this way 
one could work out all the details of the high-/c spectra P(k, z) from dj(£), once the linear 
part of P(k, z) is known. However, the need to add information on baryonic and dynamical 
parameters is not much welcome. 

Here we aim to test whether we can do better if up to 28 indipendent equations are 
allowed, a priori. In particular, no input concerning baryonic and/or dynamical parameters 
will be sought, any more. 




(5.5) 



A 



6. The SVD technique 



The risk of singularity for the matrix 

Mat = w r W i(yA) {x r )W j ( A) (x r )u r /[(3x r R(x r )] (6.1) 

arises from the vanishing of Wi - for low i values - when u r -x r attain values where Wi - for 
high i values - is still significant. Let us then introduce the SVD technique which, first of 
all, allows us to evaluate the degree of singularity of a matrix. 

The SVD technique is based on a theorem of linear algebra, stating that any real 
N r <%> N c matrix Ai, with N r > N c , can be decomposed, in a unique way apart trivial 
overall factor shifts, into a rows x columns product 

M Nr ,N c = U Nr>Nc x | diag{si) \ Nc x V^ cjA r c (6.2) 

Here, in each index site, we set the range allowed to the index there. Both Ai and V (V T 
is its transposed) are orthonormal matrices; also | diag(si) \ N is a (fully diagonal) matrix. 
If N r = N c , the inverse of Ai, in general, reads 

M~ l = V x | diag(l/si) | x U T , (6.3) 

and, as we shall discuss shortly, the technique is also a valid numerical way, not only 
to invert large matrices, but also to handle them when their full inversion is numerically 



impossible. A fair discussion of the SVD technique can be found in Numerical Recipes [18 1 



or in Matrix Computations p9|. 

The degree of singularity, however, can be inspected by considering the Sj components. 
When some of them vanishes, the matrix is singular and the problem is said to be ill- 
conditioned. Even if it is not so, as in our cases, it may happen that the ratio between 
the greatest and smallest Sj (condition number) exceeds ~ 10 6 (10 12 ), and then there is no 
hope to invert Ai in single (double) precision. Anyhow, if a level of precision 0(1 : 10 5 -10 6 ) 
is to be achieved, double precision is surely needed, so to keep the top Si/sj ratio within 
~ 10 6 -10 7 . This being necessary condition, it turns out not to be a sufficient one and, 
quite often, we need to be more restrictive. 

6.1 Singularity level 

With 28 equations available our first attempt was to invert a 28 2 matrix. As expected, its 
degree of singularity depends on the choice of u and j3. Best results are obtainable, as we 
discuss in Appendix A, when the top z r value is ~ 1.40 and the low-z domain is sampled 
above z ~ 0.05. 

In this paper, we shall consider the case u = 351, /3 = 1.8424. In Figures |6| and [/] 
we plot the Gaussian weights vs. the related redshift values and the components of the 
diagonal matrix s obtained. Let us remind that each w r is multiplied by e Xr , so that also 
w r ~ 10 -40 are not negligible. On the contrary, the Sj components do span more than 13 
o.o. m. and inverting the Ai matrix, as it is, would therefore require unusual computational 
resources. 
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Figure 6: Gauss-Laguerre weights vs. the related redshift values. 
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Figure 7: The elements Sj of the diagonal matrix, yielding the singularity level of the matrix to 
be inverted. 

The SVD approach however tells us how to deal with this case keeping as much infor- 
mation as possible from the matrix while reducing its effective singularity level. We can 
do so simply by setting to zero the unwanted Sj values. For instance, if we set to zero the 
last 7 (8) Si, the residual component are within 9 (8) orders of magnitude. The rational of 
this reduction of the condition number is explained in the cited manuals and discussed in 
Appendix A; it does not coincide with simply disregarding a number of equations, as the 
power spectra are actually recovered at all 28 redshifts, although discrepancies from input 
fluctuation spectra are safely small for 0.1 ~ z ~ 1.2. As a matter of fact, this is the very 
range where shear signals are reliable, and the matrix inversion technique substantially 
yields a best-fit recovery of input spectra. 

There is however another point deserving to be soon outlined. Real data on Cij(£) 
spectra will come with error bars. For practical application this is the critical point. In 
order to study their effect, we shall therefore add a random Gaussian noise to the Cij(£) 







3 -5 



-10 




spectra obtained by integrating the fluctuation spectra. In the inversion process, noise 
tends to be magnified: in Appendix B we show that an upper limit to noise magnification 
is approximately set by the ratio between the top and bottom Sj which were not set to 
zero. Therefore, using Sj spanning 8 orders of magnitude is likely to yield Cij(£) spectra 
completely covered by noise. 

This is a reason why we shall need to exploit the SVD technique at its limits, by keeping 
just a few non-zero S{. Lukely enough the transfered noise, in this case, is significantly 
below the theoretical upper limit and can be furtherly reduced by filtering. 

In the next subsections we shall provide quantitative details on these points. 

6.2 Recovery in the absence of noise 

Two diffeerent integration technique have been considered here: (i) a Riemannian tech- 
nique (RI, hereafter) allowing any wanted precision; (ii) a Gauss-Laguerre technique (GL, 
hereafter), that opens the door to the inversion. 

Although using 28 points and coefficients, GL does not fully succeed to recover the 
Cij(£) spectra obtained by RI. As a matter of fact, the spectra from simulations are not 
completely smooth and keep small fluctuations even for small k increments. Even a 28-th 
degree polinomial is unable to follow all of them. 

In Figure |8| we compare the results of GL and RI integration for C\\(l) and C22W 
(the coefficient H 4 is omitted). These plots are aimed to magnify the differences, which are 
0.1% at most, but can be clearly perceived. Notice, in particular, that the GL integrals 
exhibit tiny oscillations. When GL integration is performed by using a more usual number 
of points (~ 10), these oscellations risk to be confused with the effect of BAO's in fluctuation 
spectra. A similar danger exhists if we perform a RI integration after an interpolation of 
type 1 with not enough spectra along z. In this case BAO's are an effect of the interpolating 
polinomials and disappear when we increase the number of spectra along z . 



-2 



-1 1 

Log(k) 



t — i — | — i — i — i — i — | — i — i — i — i — | — i — i — i — i — | — i — r 




Log(k) 

Figure 9: Recovery of P(k, z) spectra, at the 5 redshifts shown in the frames (chosen for the sake 
of example), by using GL shear spectra and 21 (upper frame) or 20 (lower frame) non-zero s$. Black 
curves are the input spectra. The red curves are the spectra recovered; in some cases, namely for 
20-Si, the recovery is so efficient that the red lines are almost canceled by the black input spectra. 

In Figure ^ we show the results of matrix inversion if applied to GL spectra, using 20 
or 21 Si values. The passage from 9 to 8 o.o.m. in the Sj selection appears to be critical. 
Notice also that the algorithm has however some difficulty to recover the BAO shape, 
namely at low z. In the 21-Sj case, we see the low— z spectra covered by the numerical 
noise. 

In Figure [9| and in the following ones we plot a selection of spectra, for the 2j values 
shown in the frame. For the sake of homogeneity, redshifts are the same in all Figures. 
Notice that, rather than selecting those redshifts where the inversion algorithm performs 
best, we provide equi-spaced values starting from a range where inversion just begins to 
work. Among the values plotted, there are ~ 20 (mostly unplotted) redshift values where 
spectra are recovered in a fashion similar to the best results shown. 

Let us now test the results of applying the inversion algorithm, based on GL integration, 
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Figure 10: Recovery of P(k, z) spectra, at the 5 redshifts shown in the frames (the same as in the 
previous Figure), by using RI shear spectra and 16 (upper frame) or 15 (lower frame) non-zero Sj. 
Black and red curves as in the previous Figure. 



on RI spectra. In this case, the transition from unsuitable to fair spectra occurs at the 
passage from 16 to 15 values. Unsatisfatory inversion is however appreciable here through 
curve confusion, rather than by curves obscured by numerical noise. The inversion is 
slightly less precise, typical errors being 0(1 %), apart of the BAO range, where recovery 
still appears more critical; the degree of precision will be debated in detail when in the 
presence of noise. 

Before concluding the discussion on inversion in the absence of noise, let us consider 



also the case of keeping just 6 Sj, shown in Figure 11 



This Figure witnesses the high efficiency of the SVD technique, allowing to recover 
a large deal of spectra by apparently using just 6 equations. As a matter of fact, in the 
redshift range 0.2-1, fluctuation spectra are recovered with errors keeping 0(1%) in most 
of the k range; unfortunately, however, BAO's tend to be badly recovered and the low-fc 
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Figure 11: Recovery of P(k, z) spectra, at the same 5 redshifts of previous Figures, starting from 
RI shear spectra and keeping just 6 non-zero s\. Black and red curves as in the previous Figures. 

range is significantly problematic. 

6.3 Recovery in the presence of noise 

The study of error propagation from (£) to P(k, z) spectra is one of the main contribution 
of this work. As above mentioned, in the process of matrix inversion errors expand. We 
shall however show that, suitably modeling available parameters, this expansion can be 
kept under control. The case considered here should be however considered just as a good 
example. Further improvement is likely to be still possible, by refining the u-f3 selection or 
passing to a more sophisticated mapping from u and x variables. Also the noise initially 
found in recovered P(k, z) spectra will be filtered here in a simple fashion. More efficient 
filtering will be certainly possible. 

Moreover, while the 7-bin case is already significant, the huge improvement found by 
passing from 5 to 7 bins suggests a further increase in the bin number, at the limits allowed 
by shot noise. 

Bearing all that in mind, let us consider Figure [L^. Here we started from Cij(£) spectra 
obtained from an "exact" Riemann integration. From it we derived the "noisy" spectrum 

ciP(£) = C lj (£)x[l + eG}. (6.4) 

Here G is a random variable distributed in a Gaussian way around zero, with unit variance; 



e then sets the effective m.s.d.; in Figure 12 results obtained for e = 0.8% will be shown 



(cyan curces). Results are then filtered by using a 10-point top-hat filter (red curves). 
The basic SVD inversion is achieved here by keeping just 6 non-zero Sj. This allows us to 
keep C(~ 300) the ratio between top and bottom non-zero Si (condition number) and this 
is the theoretical upper limit to error magnification. In face of these predictions, we find 
that, in the central z-range and after filtering, errors can be reduced to ~ 5%, therefore 
keeping the typical error magnification well within a factor ~ 10 . This unexpected success 
is one of the results of this analysis. 
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Figure 12: Recovery of P(k, z) spectra, at the same 5 redshifts of previous Figures, starting from 
RI shear spectra and keeping just 6 non-zero Sj. The black curves, as in previous Figures, show 
the input fluctuation spectra. The rsults of inverting noisy shear spectra are the cyan curves. By 
applying to them a simpleminded top-hat filter, averaging among 10 nearby k values, errors are 
smoothed and we obtain the red curves. 



In the rest of this Section we show the errors by plotting the ratios between recovered 
and input P(k, z) for 12 equispaced values of z. We shall consider the SVD inversions 
obtained by using 7, 6, and 5 non-zero s,. The Figures [l3| these ratios are shown. In each 
Figure, 3 redshift values (shown aside) are considered by using 3 different colors. Each plot 
can be associated to the related redshift by reminding that, when z increases, the /c-range, 
for which P(k, z) is recovered, shifts to smaller values. The recovered spectra where again 
smoothed by averaging over 10 nearby points. This rough filter already yields a significant 
error reduction. More sophisticated filters could certainly allow to gain an extra factor 2, 
at least. 

As expected, reducing the number of non-zero Sj reduces error magnification. In order 
to achieve a tolerable error magnification, however, the number of non-zero s, needs to 
be small. According to Figure |f| taking 6 non-zero (over 28) appear be the best 
compromise: further reducing the non-zero Sj number starts to affect significantly P(k, z) 
normalization; keeping more non-zero Sj furtherly magnifies errors. 

This is also confirmed by Figure 14, which shows the average shifts and the variance 
about them for the points 2-27. This Figure shows a sort of oscillatory residuals, reminis- 
cent of the use of a finite number of orthogonal polynomials. Let us also notice that, if 
an improved filtering furtherly reduces errorbars, using 7 non-zero Sj could also be a fair 
choice, as oscillatory residuals decrease when the number of non-zero Sj is greater. 

With such low numbers of non-zero Sj, recovering the BAO shape appears however 
difficult; in the low /c-range the inversion seems to work poorly. 



7. Discussion 



In this paper we discuss the invertion of the integral expression yielding the tomographic 
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Figure 13: Ratio between recovered and input power spectra at different redshifts (at the r.h.s. of 
the Figures), when the shear spectra is added a Gaussian noise ~ .8 %. The recovered spectra where 
filtered through a 10-point top-hat smoothing. The Figures are obtained when 7, 6, 5 non-zero Si 
are kept. Reducing the number of non-zero s, reduces error magnification, but biases the inversion 
process. In the case considered, best results seem obtainable with 6 non-zero s$. See also the next 
Figure. 
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shear spectra C{j{£) from the density fluctuation spectra P(k, z). In a previous work (Paper 
I), the question had been discussed with (N=)5 bin tomography, yielding N(N+1)/2=15 
independent shear spectra. With present samples, spreading shear data over 5 bins yields 
spectra already seriously affected by shot noise. We however expect future data (Euclid) 
to cover a larger sky portion, so allowing up to 10 bin tomography. 

We however kept below this limit and considered here N=7, yielding N(N+l)/2=28 
independent shear spectra. This already allows us great improvements in respect to N=5. 
Another contribution of this work is the study of noise propagation from tomographic shear 
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Figure 14: Average displacements from unity of the ratios between input and recovered spectra; 
m.s. oscillation amplitudes around average are shown by errorbars. 

data to fluctuation spectra. 

We have started from the equation yielding shear from fluctuation spectra in the form 



c A = Ma,tPt 



(7.1) 



Here Cij{l) = H 4 ca(£) and p r (i) = P(u r /£,u r ), u being the comoving distance of the 
redshift z. Let us remind that fluctuation spectra used here were obtained by suit- 
ably connecting linear and hydro simulation spectra, according to a technique widely 
discussed in Paper I. The basic issue here is that we need no approximate spectral ex- 
pression (like HALOFIT) to perform a fairly normalized connection, thanks to the wide 



box (410/i _1 Mpc) used in simulations. Eq. (7T) deals with each I separately; summation 
replaces there integration and is performed by using a Gauss-Laguerre (GL) technique, 
involving the replacement of the comoving distance variable u with 



(u/uY 



(7.2) 



and approximating the large-u behavior of integrand functions by exp(— x). As A takes 28 



values, we are allowed a 28-point GL integration; the sum (7.1) then closely approaches the 
exact integral. Tiny residual differences however remain, which will matter in the discussion 
below. The matrix Ai is obtained from the expected galaxy distribution on z, the bin 
structure, the photometric-spectrographic redshift relation, and the space geometry. 

Recovering P(u r /£,u r ) from Cij(£) then requires inverting the matrix Ai. Unfortu- 
nately, as different bins refer to different depths, Ai approaches a singular behavior. 

The recovery of P(k, z) is therefore based on the singular value decomposition (SVD) 
technique for matrix inversion. Such technique soon yields the degree of singularity through 
the elements of a 28<%>28 diagonal matrix s. Singularity however depends on the choice of 
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the parameters u and (3 and, in this work, we provide an analysis of the criteria for their 
selection. No full Ai matrix can be however inverted with standard numerical means, as 
the Si diagonal components of the matrix s span ~ 13-14 o.o.m., in spite of optimizing 
the u-(3 selection. (A look at the Figure in Appendix A however shows that an unsuitable 
parameter choice can yield s; spanning many more o.o.m. 's.) 

According to the SVD technique, we can however reduce the divergence level, by setting 
to zero a number of Sj. This is similar - but not fully equivalent and, indeed, not so drastic - 
to suppressing a number of equations. Our analysis greatly exploits this option, also in 
view of the fact that also error propagation depends on the ratio between minimum and 
maximum non-zero Sj. 

The simplest task that the inversion procedure can accomplish is the recovery of P(k, z) 
from Cij(£) obtainable via GL integration. The only limitation that we meet, in this case, 
is the divergence degree outlined by the diagonal s matrix components. Out of 28 Sj 
components, we find that best results are obtainable when keeping 20 non-zero Sj. Already 
with 21 components, numerical noise affects the low-z recovered spectral components. 

We then tested the inversion of the results of exact integration. Although being closer 
to realistic, this works as though adding a noise to the ca(£) components. As a matter 
of fact, to obtain fair results, we must then furtherly reduce the number of non-zero Sj. 
Best results are obtained with 15 Sj components. Typical discrepancies between input and 
output fluctuation spectra are then C(l %), apart of some wider discrepancy arising in the 
low-A; BAO range, namely for low redshift. 

These results are however a great improvement in respect to 5 bins. Then, the very 
normalization of recovered spectra was at risk and fair results could be obtained only 
through renormalization in the low-A: range. 

We then tested error propagation by superimposing a Gaussian noise (m.s.a. ~ 0.8%) 
onto the Cij(£) spectra obtained from the integration of simulation spectra. Such noise is 
however magnified by matrix inversion. A theoretical upper limit on noise magnification 
is set by the ratio between the top and bottom Sj kept (see Appendix A). This calls for 
pushing further down the number of non-zero Sj. As a matter of fact, best results are 
obtained when keeping 6 non-zero Sj (over 28!). If 7 are kept, the error is greater. If we 
try to keep just 5 of them, we meet displacements on recovered spectrum normalization. 

The noise in direct inversion results was then tentatively filtered by using a 10-point 
top-hat smoothing. After using this admittedly rough smoothing technique, we find that 
the overall noise magnification is less than a factor 10, in respect of a theoretical upper limit 
~ 300 . Let us also outline that spectral values on nearby k values derive from independent 
i inputs, so that filtering acts in a "normal" direction, in respect to matrix inversion. 

The /c-region most penalized by the suppression of a large number of Sj is the BAO 
range. Here errors are greater and we can summarize our outputs by saying that this 
technique does not allow to recover the BAO structure. 

8. Conclusions 

Shear data are expected to be a new and effective resource to discriminate among cosmo- 



logical models. They are most sensitive to the rate of expansion and fluctuation growth 
at low— z, when the contribution of Dark Energy is determinant. Accordingly, they are 
expected to be an outstanding probe on the DE equation of state w{z) and/or to possible 
deviation from GR, or to couplings between the dark cosmic components. 

From data we expect tomographic shear spectra to be derived. A basic difficulty 
amounts to cleansing rough data from intrinsic ellipticities. This fundamental question is 
relevant for our present aims just because it might constrain the number N of tomographic 
bins used. In fact, to our aims, using TV = 7 is found to yield much better results than 
N = 5 or less. 

An analysis of spectral data can be performed in accordance with a Bayesian paradigm. 
In this context, Cij(£) model predictions will be made, and compared with Cij(£) data and 
their errorbars. 

Other sets of data, possibly derived from the same Euclid experiment, will directly 
concern the P(k, z) spectra. 

The possibility to skip Cij(£) model predictions, by directly using the P(k,z) spectra 
obtained from them, together with other data on P(k, z), could then be appealing, although 
it is hard to state, at this stage, if error magnification is an obstacle to this option. 

In order to build the A4 matrix, that shall be inverted to work out P(k, z) from 
Cij(£), we must be able to convert redshifts into comoving distances. This requires some 
information on the model geometry, but no dynamical or baryonic information: neither 
as and n s , nor parameters describing baryon physics (including f^) or cosmic reionization 
shall be input. 

For the sake of example, let us then outline a possible investigation pattern: SNIa data 
could become so good that the distance modulus could provide u{z) with negligible errors. 
From such u(z) the parameters of the background cosmology can be derived. However, 
independently from such derivation, we could directly use u(z) to build the window func- 
tions W r and derive from them the inverted matrix M. -1 , so obtaining P(k,z) -apart of 
a constant factor (VLmH^) 2 - and the fluctuation growth law G(a). In parallel with fitting 
u(z) to models, we then also fit G{a) to them, so performing two independent model tests, 
based on geometry and dynamics, respectively. Altogether, the two tests could provide 
a direct confirm/falsification of models where DE and DM are two real and independent 
cosmic components. 

More in general, this technique allows us to measure the fluctuation growth in an unbi- 
ased way, over any fc-range, namely those where different hypotheses on stellar formation, 
SN explosion, AGN energy release, etc., are critical. Accordingly, the inversion technique 
might become an important tool to discriminate among different hypotheses on baryon 
physics, possibly testing parameters beyond those spanning the bayesian parameter space. 

Let us finally outline that there is a substantial space for further improvements of the 
technique providing P(k, z) from shear spectra. We wish to outline, in particular, the fol- 
lowing points: (i) We saw the great improvement achieved when passing from N=5 to N=7; 
here we did not test greater N values, but it is licit to expect that further improvements 
are at hand if using 8, 9, or 10 bins is possible. Although Euclid is expected to provide 
a tomography with up to 10 bins, it is not clear to us up to which point intrinsic shear 



cleansing is more effective with N = 5 or less, (ii) Noise propagation was considered here 
to test error propagation. Propagated errors were then filtered, with a procedure effective 
but rought. Better filtering techniques are surely available, but selecting among them de- 
pends on the actual nature of the "noise" to be considered, (iii) A specific technical issue 
concerns the passage from the the comoving distance u to the integration coordinate x. 
Changes in the parameters u and ft entering in the conversion expression x = (u/u) 13 were 
found to be critical, and an optimal parameter choice to be almost determinant. Let us 
then outline that u—x mappings more effective than a simple power low, namely in the 
low-x range, could also allow us further substantial improvements. 
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Appendix A 

1. The singular value decomposition (SVD) 

The SVD technique is based on a theorem of linear algebra stating that any M x N 
(M > N) matrix M can be decomposed in the matrix product U ■ s ■ V T (here V T is the 
transpose of a matrix V) or, more in detail, 

Mu ... M N i 

Mim ■■■ Mnm 
with U and V being orthogonal, i.e., 

UiaUib = Sab , ^ V iaVib = &ab , {A2) 

i i 



Umi 



U N1 



Unm 
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v NN 
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or 

U T ■ U = 1 , V T ■ V = 1 (A3) 

and s being diagonal. 

In this work we are interested just in the case when M = N and all matrices are 
square; accordingly, it is also U ■ U T = 1 (and V ■ V T = 1). 

The decomposition (Al) is unique, apart trivial multiplicative factors. 

The condition number of the matrix M is then defined as the ratio K between the 
largest and smallest (in magnitude) Sj. A matrix is said to be ill-conditioned if K is too 
large (or singular if K — > oo). 

2. Norms and error magnification 

Let us remind first that the (2-)norm of a vector x reads 

1/2 



2Ci I 



Xi being its components. The (2-)norm of a matrix A is then defined as follows: 

|A| = max (A4) 

i.e.: | A | is the maximum factor by which the matrix A can amplify a non-vanishing vector 
x. It can then be shown that, if we perform a SVD of A and s\ is the top component of 
the diagonal matrix s, in the relation A = U • s • V" 1 , it is 

|A| = si . 

Let us then outline that 

A" 1 = V-t-U- 1 , 

with ti = s^ 1 . Accordingly, \A~ 1 \ = l/s n , s n being the smallest component of the matrix 
s and 

|A| • |A~ X | = Sl /s n . (AS) 
The matrix-vector and matrix-matrix products have then the following properties: 

| A • x| < |A| • |x| (A6) 

|A-B| < |A| • |B| . 

Let us then assume that 

A • x = b (A7) 

and consider x = x + Jx, as we do when we try to reobtain P(k, z) from noisy Cij(£). It 
shall be A • x = b + 5h and therefore 



so that, owing to eq. (A6), 
In turn, eq. (A7) yields 

and, therefore 



|<5x| < • \5b 



1 |A| 
W ~ TbT 



Idx| < lAI • IA- 1 ! I5b| 



|x| ~ 1 1 1 1 |b| 
Owing to eq. (A5) we have then that 



|(5x| si \5b\ , ._. 

— TuT • 



x 



i.e., that error magnification is surely smaller than the ratio between the maximum and 
minimum s components considered, coinciding with the condition number K. 



For more details on this Section and proofs of the theorems, besides of [18, 16], see 
also fl2C 



3. Gauss Laguerre parameter selection 

The Gauss-Laguerre sum, replacing the integration over P(k, z) to obtain the shear 
spectra CL- (£) , ought to by made with an eye to the matrix inversion which will follow. The 
condition number K of the matrix M, in fact, depends in a critical way on the selection 
of the u and (3 parameters used in the change of variable x = (u/u)? . A bad selection of 
u and (3 can yield K values up to 10 20 or even larger. Moreover, K is also directly linked 
to error magnification, according to eq. (A7). 

To minimize K we can act on three parameters: besides of (3 and u, we must take into 
account N cu t, the number of Si which (unavoidably) will be set to zero. Increasing N cu t 
obviously reduces K, but our aim amounts to minimize K while keeping N cu t sufficiently 
low, so to grant a fair and physical solutions to the a substantial part of the 28 equations 
forming the linear system. 

Let us then keep into account that the galaxy redshift distribution n(z) rapidly declines 
beyond z ~ 1 and that, therefore, the window functions W r and the M matrix element 
also fade. In the following we shall not directly consider contributions coming from above 
a maximum redshift z max = 1.4, corresponding to a maximum comoving distance u ma x = 
4194 Mpc, a choice yielding fair results. We shall then select the value of u so that 



MP 



xm being the largest node in the Gauss-Laguerre integration. 

With this rule, we can then estimate the condition number for a wide range of f3 (from 
0.1 to 3.5) and N cu t values. The latter parameter will be considered in the interval (we 
keep all singular values) - 26 (keeping only the two largest singular values to estimate their 
ratio). 



Figure 15: The condition number K is plotted against /3 and the number of non-zero Si kept. 
At high K, there appear some graphic irregularities in the curves, due to the discreteness of the 
number. 



In Figure 15, K is plotted as a function of (3 and of the residual Si number (28 — N cut ). 
The plot shows how dramatically K blows up just by choosing unsuitable values of the free 
parameters. It is also evident that, reducing the number of non-zero Sj, K sistematically 
decreases and that, for a given number of non-zero Si, K depends on f3. 

However, when the number of non-zero Si is small (<~ 10) the K dependence on f3 
weakens. This can be exploited to perform slight changes of the ft value, yielding a suitably 
different set of u r for the Gauss-Laguerre nodes x r . In this way, we can obtain P(k) at 
almost any z value, in the z-range where matrix inversion is efficient (see text). 



